#ifndef SPHERE_H
#define SPHERE_H

#include <vector>
#include <cmath>
#include <iostream>
#include <string>

// 球极投影：从球面到平面（单个double参数）
void stereographicProjection(double x, double y, double z, double& X, double& Y) {
    if (z == 1) {
        X = Y = NAN;  // 北极点特殊处理
    } else {
        X = x / (1 - z);
        Y = y / (1 - z);
    }
}

// 球极投影：从球面到平面（处理向量）
void stereographicProjection(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, 
                              std::vector<double>& X, std::vector<double>& Y) {
    size_t n = x.size(); // 假设 x, y, z 的大小相同
    X.resize(n);
    Y.resize(n);

    for (size_t i = 0; i < n; ++i) {
        if (z[i] == 1) {
            X[i] = Y[i] = NAN;  // 北极点特殊处理
        } else {
            X[i] = x[i] / (1 - z[i]);
            Y[i] = y[i] / (1 - z[i]);
        }
    }
}

// 反向投影：从平面到球面（单个double参数）
void inverseStereographicProjection(double X, double Y, double& x, double& y, double& z) {
    double r = std::sqrt(X * X + Y * Y);
    z = (r * r - 1) / (r * r + 1);
    x = 2 * X / (r * r + 1);
    y = 2 * Y / (r * r + 1);
}

// 反向投影：从平面到球面（处理向量）
void inverseStereographicProjection(const std::vector<double>& X, const std::vector<double>& Y, 
                                    std::vector<double>& x, std::vector<double>& y, std::vector<double>& z) {
    size_t n = X.size();
    x.resize(n);
    y.resize(n);
    z.resize(n);

    for (size_t i = 0; i < n; ++i) {
        double r = std::sqrt(X[i] * X[i] + Y[i] * Y[i]); 
        z[i] = (r * r - 1) / (r * r + 1);
        x[i] = 2 * X[i] / (r * r + 1);
        y[i] = 2 * Y[i] / (r * r + 1);
    }
}

// 计算两点间的欧几里得距离
double euclideanDistance(double x1, double y1, double x2, double y2) {
    return std::sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
}

// 计算累积弦长并归一化
void calculateCumulativeChordLengths(const std::vector<double>& x, const std::vector<double>& y,
                                     std::vector<double>& t, std::vector<double>& x_t, std::vector<double>& y_t) {
    int N = x.size();  // 数据点的数量

    // 初始化
    t.resize(N);       // 存储累积弦长
    x_t.resize(N);     // 存储对应的 x(t_i)
    y_t.resize(N);     // 存储对应的 y(t_i)

    // 第一个点的累积弦长是 0
    t[0] = 0.0;
    x_t[0] = x[0];
    y_t[0] = y[0];

    // 计算其他点的累积弦长
    for (int i = 1; i < N; ++i) {
        // 计算当前点与前一个点的欧几里得距离
        double dist = euclideanDistance(x[i-1], y[i-1], x[i], y[i]);

        // 累积弦长
        t[i] = t[i-1] + dist;

        // 存储 x(t_i) 和 y(t_i)
        x_t[i] = x[i];
        y_t[i] = y[i];
    }

    // 找到累积弦长的最大值
    double max_t = t[N-1];

    // 归一化累积弦长 t，使其范围为 [0, 1]
    for (int i = 0; i < N; ++i) {
        t[i] /= max_t;  // 将 t[i] 归一化
    }
}

// 计算任意维度空间中的欧几里得距离
double euclideanDistance(const std::vector<double>& p1, const std::vector<double>& p2) {
    // 假设 p1 和 p2 的维度相同
    double dist = 0.0;
    for (size_t i = 0; i < p1.size(); ++i) {
        dist += (p2[i] - p1[i]) * (p2[i] - p1[i]);
    }
    return std::sqrt(dist);
}

// 计算累积弦长并归一化
void calculateCumulativeChordLengths(const std::vector<std::vector<double>>& points,
                                     std::vector<double>& t, 
                                     std::vector<std::vector<double>>& points_t) {
    int N = points.size();  // 数据点的数量

    // 初始化
    t.resize(N);             // 存储累积弦长
    points_t.resize(N);      // 存储对应的点(t_i)

    // 第一个点的累积弦长是 0
    t[0] = 0.0;
    points_t[0] = points[0];

    // 计算其他点的累积弦长
    for (int i = 1; i < N; ++i) {
        // 计算当前点与前一个点的欧几里得距离
        double dist = euclideanDistance(points[i-1], points[i]);

        // 累积弦长
        t[i] = t[i-1] + dist;

        // 存储点
        points_t[i] = points[i];
    }

    // 找到累积弦长的最大值
    double max_t = t[N-1];

    // 归一化累积弦长 t，使其范围为 [0, 1]
    for (int i = 0; i < N; ++i) {
        t[i] /= max_t;  // 将 t[i] 归一化
        //std::cout << "t[i]=" << t[i] << std::endl;
    }
}

#endif